#include <iostream>
#include <string>
#include <vector>
#include <complex>
#include <fstream>
#include <iterator>
#include <numeric>
#include <algorithm>

using namespace std;

struct power2
{
    double operator()(const complex<double>& arg)
    {
        double a = abs(arg);
        return a * a;
    }
};
                
class GaussReader
{
public:
    GaussReader(std::string fileName)
        :
            inp(*new ifstream(fileName.c_str())) {}
    
    double getValue()
    {
        double rslt;
        inp >> rslt;
        return rslt; 
    }
    
    ~GaussReader()
    {
        delete &inp;
    }
    
private:
    ifstream& inp;
};

template
<
    typename SignalSystem
>
double calcEnergy(const SignalSystem& signals)
{
    double energy = 0.0;
    
    for (unsigned i = 0; i < signals.size(); ++ i)
    {
        for (unsigned j = 0; j < signals[i].size(); ++ j)
        {
            double c = abs(signals[i][j]);
            energy += c * c;
        }
    }
    
    return energy;
}

template
<
    typename SignalSystem
>
double calcBer(SignalSystem& signals, double snr, unsigned count = 1000)
{
    GaussReader gauss("gauss.txt");
    
    double energy = calcEnergy(signals);
    
    double errors = 0.0;
    
    for (unsigned q = 0; q < count; ++ q)
    {
        for (unsigned i = 0; i < signals.size(); ++ i)
        {
            auto signal = signals[i];
            
            for (unsigned j = 0; j < signal.size(); ++ j)
            {
                double x = pow(10.0, - snr / 20) * energy * gauss.getValue();
                double y = pow(10.0, - snr / 20) * energy * gauss.getValue();
                signal[j] += complex<double>(x, y);
            }
            
            double sumopt = 0;
            unsigned jopt = 0;
            
            for(unsigned j = 0; j < signals.size(); ++ j)
            {
                auto diff = signal;
                
                transform(
                    signal.begin(), signal.end(),
                    signals[j].begin(), diff.begin(),
                    minus<complex<double>>()
                );
                
                vector<double> rslt(diff.size());
                
                transform(
                    diff.begin(), diff.end(),
                    rslt.begin(),
                    power2()
                );
                
                double sum = accumulate(
                    rslt.begin(), rslt.end(),
                    0.0
                );
                
                if (j == 0)
                    sumopt = sum;
                else
                    sumopt = min(sum, sumopt);
                
                if(sumopt == sum)
                    jopt = j;
            }

            if (jopt != i)
                ++ errors;
        }
    }
    
    return errors / (count * signals.size());
}

int main(int argc, char** argv)
{
    using namespace std;
 
    cout << pow(10, 1) << endl;   
    typedef vector<complex<double>> Signal;
    typedef vector<Signal> SignalSystem;
    
    SignalSystem bpsk =
    {
        { complex<double>(1, 0) },
        { complex<double>(-1, 0) }
    };
    
    cout << bpsk[0][0] << endl;
    cout << bpsk[1][0] << endl;
    
    vector<double> snr = {-20, -10, -5, 0, 5, 10, 20};
    vector<double> result;
    
    for (auto it = snr.begin();
        it != snr.end();
        ++ it)
    {
        double ber = calcBer(bpsk, *it, 1000000);
        result.push_back(ber);
        cout << ber << endl;
    };
}
